#GOALS
#1.Describe how to built UPSS.
#2.Exploring UPSS effects on canopy microclimate, including temperature and light intensity, relative to the open and shrub.
#3.Understanding how different light permeabilities and shelter shapes influence the above parameters.
library (tidyverse)
## Warning: package 'tidyverse' was built under R version 3.6.2
## -- Attaching packages ----------------------- tidyverse 1.3.0 --
## v ggplot2 3.2.1 v purrr 0.3.3
## v tibble 2.1.3 v dplyr 0.8.3
## v tidyr 1.0.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.4.0
## Warning: package 'forcats' was built under R version 3.6.2
## -- Conflicts -------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(ggplot2)
#import datasets
panoche.climate.2019<-read.csv("C:/Users/Nargol Ghazian/Desktop/Animal-Behaviour-and-Climate-project/CH3/panoche_climate_hourly 2019.csv")
mean.daily.panoche <- read.csv("C:/Users/Nargol Ghazian/Desktop/Animal-Behaviour-and-Climate-project/CH3/weather station climate.csv")
str(panoche.climate.2019)
## 'data.frame': 8760 obs. of 28 variables:
## $ station.id : int 56 56 56 56 56 56 56 56 56 56 ...
## $ site : Factor w/ 2 levels "","panoche": 2 2 2 2 2 2 2 2 2 2 ...
## $ station.name: Factor w/ 2 levels "","Los Banos": 2 2 2 2 2 2 2 2 2 2 ...
## $ CIMIS.Region: Factor w/ 2 levels "","San Joaquin Valley": 2 2 2 2 2 2 2 2 2 2 ...
## $ hour : int 1 2 3 4 5 6 7 8 9 10 ...
## $ date : Factor w/ 27 levels "","5/18/2019",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ month : int 5 5 5 5 5 5 5 5 5 5 ...
## $ day : int 18 18 18 18 18 18 18 18 18 18 ...
## $ year : int 2019 2019 2019 2019 2019 2019 2019 2019 2019 2019 ...
## $ precip : num 0 0 0 0 0 0 0 0 0 0 ...
## $ radiation : int NA NA NA NA NA 57 176 432 597 723 ...
## $ air.temp : num 8.7 9.2 8.8 7.6 6.6 7 10.3 13.7 16 17.5 ...
## $ wind.speed : num 1.4 1.3 0.9 1.2 1.2 0.6 1 0.9 1.4 3 ...
## $ soil.temp : num 18.1 18 18 17.9 17.8 17.7 17.6 17.5 17.4 17.3 ...
## $ air.temp.F : num 47.7 48.6 47.8 45.7 43.9 ...
## $ X : logi NA NA NA NA NA NA ...
## $ X.1 : logi NA NA NA NA NA NA ...
## $ X.2 : logi NA NA NA NA NA NA ...
## $ X.3 : logi NA NA NA NA NA NA ...
## $ X.4 : logi NA NA NA NA NA NA ...
## $ X.5 : logi NA NA NA NA NA NA ...
## $ X.6 : logi NA NA NA NA NA NA ...
## $ X.7 : logi NA NA NA NA NA NA ...
## $ X.8 : logi NA NA NA NA NA NA ...
## $ X.9 : logi NA NA NA NA NA NA ...
## $ X.10 : logi NA NA NA NA NA NA ...
## $ X.11 : logi NA NA NA NA NA NA ...
## $ X.12 : Factor w/ 3 levels " ","R","Y": 1 1 1 1 1 1 1 1 1 1 ...
na.omit (panoche.climate.2019)#exclude missing values
## [1] station.id site station.name CIMIS.Region hour
## [6] date month day year precip
## [11] radiation air.temp wind.speed soil.temp air.temp.F
## [16] X X.1 X.2 X.3 X.4
## [21] X.5 X.6 X.7 X.8 X.9
## [26] X.10 X.11 X.12
## <0 rows> (or 0-length row.names)
head(is.na(panoche.climate.2019))#check for missing values
## station.id site station.name CIMIS.Region hour date month day
## [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [4,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [5,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [6,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## year precip radiation air.temp wind.speed soil.temp air.temp.F X
## [1,] FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE
## [2,] FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE
## [3,] FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE
## [4,] FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE
## [5,] FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE
## [6,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
## X.1 X.2 X.3 X.4 X.5 X.6 X.7 X.8 X.9 X.10 X.11 X.12
## [1,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [2,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [3,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [4,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [5,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [6,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
head(na.omit (mean.daily.panoche))
## station.name site sensor lat long date year month day
## 1 los banos panoche station -120.867 37.0563 5/18/2019 2019 5 18
## 2 los banos panoche station -120.867 37.0563 5/19/2019 2019 5 19
## 3 los banos panoche station -120.867 37.0563 5/20/2019 2019 5 20
## 4 los banos panoche station -120.867 37.0563 5/21/2019 2019 5 21
## 5 los banos panoche station -120.867 37.0563 5/22/2019 2019 5 22
## 6 los banos panoche station -120.867 37.0563 5/23/2019 2019 5 23
## temp.max temp.min temp.mean
## 1 83.1 55.2 69.15
## 2 83.3 55.4 69.35
## 3 83.5 55.4 69.45
## 4 83.7 55.6 69.65
## 5 83.9 55.8 69.85
## 6 84.2 55.9 70.05
#DATA VIZ MACRO (WEATHER STATION)
##macro-climate plots
ggplot (panoche.climate.2019, aes((hour), air.temp)) + geom_smooth() + xlab("Hour (0-24)") + ylab ("Temperature (°C)")+ theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))#how air temperature changed during 24h period in 2019 during study period
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 8136 rows containing non-finite values (stat_smooth).

ggplot (panoche.climate.2019, aes((day), air.temp)) + geom_smooth() + xlab("Day") + ylab ("Temperature (°C)")+ theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))#air temperature over the days of the study period
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 8136 rows containing non-finite values (stat_smooth).

ggplot (panoche.climate.2019, aes((hour), soil.temp)) + geom_smooth() + xlab("Hour (0-24)") + ylab ("Temperature (°C)")+ theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))#how soil temperature changed during 24h period in 2019 during study period
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 8136 rows containing non-finite values (stat_smooth).

ggplot (panoche.climate.2019, aes((day), soil.temp)) + geom_smooth() + xlab("Day") + ylab ("Temperature (°C)")+ theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))#how soil temperature changed during 24h period in 2019 during study period
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 8136 rows containing non-finite values (stat_smooth).

ggplot (panoche.climate.2019, aes((hour), radiation)) + geom_smooth() + xlab("Hour (0-24)") + ylab ("Solar Radiation (W/m²)")+ theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))#how solar radiation changed during 24h period in 2019 during study period
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 8347 rows containing non-finite values (stat_smooth).

ggplot (panoche.climate.2019, aes((day), radiation)) + geom_smooth() + xlab("Day") + ylab ("Solar Radiation (W/m²)")+ theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))#sunlight experinced over the the study period
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 8347 rows containing non-finite values (stat_smooth).

library (ggbeeswarm)
## Warning: package 'ggbeeswarm' was built under R version 3.6.3
ggplot (panoche.climate.2019, aes(as.factor(date), air.temp)) + geom_boxplot() + xlab("Date") + ylab ("Temperature (°C)")+ geom_quasirandom(alpha=0.05)+ theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))+ theme(axis.text.x = element_text(angle = 90))+ geom_smooth(se=FALSE, color="black", aes(group=1))#daily air temperature averages
## Warning: Removed 8136 rows containing non-finite values (stat_boxplot).
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 8136 rows containing non-finite values (stat_smooth).
## Warning: Removed 8136 rows containing missing values
## (position_quasirandom).

#DATA VIZ MICRO
##micro-climate plots
##let's try some plots for shelter, shrub, and open
shelter.shrub.open <- read.csv("C:/Users/Nargol Ghazian/Desktop/Animal-Behaviour-and-Climate-project/CH3/shrub_contrast_final.csv")#import dataset
ggplot(shelter.shrub.open, aes((microsite), temp, fill=microsite)) + geom_boxplot() + xlab("Microsite") + ylab ("Temperature (°F)")+ theme_classic()+ theme(axis.text=element_text(size=12))+theme(axis.text.x = element_text(angle = 90))+theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())+ labs(fill = "Microsite")+ stat_summary(fun.y=mean, colour="black", geom="point", shape=18, size=3,show_guide = FALSE)#Boxplot temperaure by microsite
## Warning: `show_guide` has been deprecated. Please use `show.legend`
## instead.

ggplot(shelter.shrub.open, aes((microsite), intensity, fill=microsite)) + geom_boxplot() + xlab("Microsite") + ylab ("Solar Radiation (lum/ft²)")+ theme_classic()+ theme(axis.text=element_text(size=12))+theme(axis.text.x = element_text(angle = 90))+theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())+ labs(fill = "Microsite")+ stat_summary(fun.y=mean, colour="black", geom="point", shape=18, size=3,show_guide = FALSE)#Boxplot light intensity by microsite
## Warning: `show_guide` has been deprecated. Please use `show.legend`
## instead.
## Warning: Removed 9330 rows containing non-finite values (stat_boxplot).
## Warning: Removed 9330 rows containing non-finite values (stat_summary).

ggplot(shelter.shrub.open, aes((day), temp, color=microsite)) + geom_smooth()+ xlab("Day") + ylab ("Temperature (°F)")+ theme_classic()+ theme(axis.text=element_text(size=12))+stat_summary(fun.y=max, geom="point", size=2, aes(shape = microsite))+ theme(legend.position = "none") + facet_grid("microsite")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

ggplot(shelter.shrub.open, aes((day), intensity, color=microsite)) + xlab("Day") + ylab ("Solar Radiation (lum/ft²)")+ theme_classic()+ theme(axis.text=element_text(size=12))+geom_smooth()+ stat_summary(fun.y=max, geom="point", size=2, aes(shape = microsite))+ theme(legend.position = "none") + facet_grid("microsite")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 9330 rows containing non-finite values (stat_smooth).
## Warning: Removed 9330 rows containing non-finite values (stat_summary).

library(ggpubr)
## Loading required package: magrittr
##
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
##
## set_names
## The following object is masked from 'package:tidyr':
##
## extract
ggline(shelter.shrub.open, x = "day", y = "temp", color = "microsite",
add = "mean_ci", shape = "microsite", xlab = "Day", ylab = "Temperature (°F)", legend.title= "Microsite", legend="right")#dotted line graph of mean daily temperatures

ggline(shelter.shrub.open, x = "day", y = "intensity", color = "microsite",
add = "mean_ci", shape = "microsite", legend.title= "Microsite", xlab = "Day", ylab = "Solar Radiation (lum/ft²)", legend="right")#dotted line graph of mean daily solar radiation
## Warning: Removed 9330 rows containing non-finite values (stat_summary).

#visualizing mean, median, and mode with histograms
ggplot(shelter.shrub.open, aes(temp, fill = microsite)) +
geom_histogram(binwidth = 5) +
scale_fill_brewer(palette = "Set1")+ labs(fill = "", x = "Temperature (°F)", y = "Frequency")+theme_classic()+ theme(axis.text=element_text(size=12))+labs(fill = "Microsite")

library(ggpubr)
ggplot(shelter.shrub.open, aes(temp, fill = microsite)) +
geom_histogram(binwidth = 5) +
scale_fill_brewer(palette = "Set1")+ labs(fill = "", x = "Temperature (°F)", y = "Frequency")+theme_classic()+ theme(axis.text=element_text(size=12))+ stat_central_tendency(aes(color = microsite), type = "mean", color="green", linetype = 2)+ stat_central_tendency(type = "median", color = "blue", linetype = 2)+ facet_grid(~microsite)+theme(legend.position = "none")

ggplot(shelter.shrub.open, aes(temp, fill = microsite)) +
geom_density() +
scale_fill_brewer(palette = "Set1")+ labs(fill = "", x = "Temperature (°F)", y = "Frequency")+theme_classic()+ theme(axis.text=element_text(size=12))+labs(fill = "Microsite")

ggplot(shelter.shrub.open, aes(temp, fill = microsite)) +
geom_density() +
scale_fill_brewer(palette = "Set1")+ labs(fill = "", x = "Temperature (°F)", y = "Frequency")+theme_classic()+ theme(axis.text=element_text(size=12))+ stat_central_tendency(aes(color = microsite), type = "mean", color="green", linetype = 2)+ stat_central_tendency(type = "median", color = "blue", linetype = 2)+ facet_grid(~microsite)+theme(legend.position = "none")

#EDA
#explore normality
#can't do Shapiro-Wilk dataset too large
library(ggpubr)
ggqqplot(shelter.shrub.open$temp)#Data are positively(right) skewed, Gaussian

hist(shelter.shrub.open$temp)

ggqqplot(shelter.shrub.open$intensity)#Data follow a poisson distribution
## Warning: Removed 9330 rows containing non-finite values (stat_qq).
## Warning: Removed 9330 rows containing non-finite values (stat_qq_line).
## Warning: Removed 9330 rows containing non-finite values (stat_qq_line).

hist(shelter.shrub.open$intensity)

#Neither intensity nor temperature follow a normal distribution
#explore light and temperature relationship
macro.micro.contrast <- read.csv("C:/Users/Nargol Ghazian/Desktop/Animal-Behaviour-and-Climate-project/CH3/micro-macro-contrast.csv")
cor.test(macro.micro.contrast$temp, macro.micro.contrast$intensity, method = "kendall")#light and temperature are correlated p<0.0001, tau=0.2813
##
## Kendall's rank correlation tau
##
## data: macro.micro.contrast$temp and macro.micro.contrast$intensity
## z = 45.612, p-value < 2.2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## 0.2659057
#scatterplot for the relationship
ggplot(macro.micro.contrast, aes(temp, intensity))+ geom_point()+ ylab("Light Intensity (lum/ft²)")+ xlab("Temperature (°F)")+geom_smooth(method='glm', se=FALSE)+theme_classic()#Perhaps too dense too look at
## Warning: Removed 9330 rows containing non-finite values (stat_smooth).
## Warning: Removed 9330 rows containing missing values (geom_point).

ggplot(macro.micro.contrast, aes(temp, intensity))+geom_point()+ ylab("Light Intensity (lum/ft²)")+ xlab("Temperature (°F)")+geom_smooth(method='glm', se=FALSE)+ facet_grid(~microsite)+ theme_classic()#faceted by microsite
## Warning: Removed 9330 rows containing non-finite values (stat_smooth).
## Warning: Removed 9330 rows containing missing values (geom_point).

ggplot(macro.micro.contrast, aes(temp, intensity))+geom_point()+ ylab("Light Intensity (lum/ft²)")+ xlab("Temperature (°F)")+geom_smooth(method='glm', se=FALSE)+ facet_grid(~cover.type)+ theme_classic()#faceted by cover type
## Warning: Removed 9330 rows containing non-finite values (stat_smooth).
## Warning: Removed 9330 rows containing missing values (geom_point).

#MODELS
##model hypotheses
##pair-wise comparison
library(sjPlot)
## Warning: package 'sjPlot' was built under R version 3.6.3
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
lm.temp <- glm(temp~as.factor(microsite), data = shelter.shrub.open, family="gaussian")
summary (lm.temp)
##
## Call:
## glm(formula = temp ~ as.factor(microsite), family = "gaussian",
## data = shelter.shrub.open)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -39.543 -16.496 -5.442 13.796 110.218
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 73.6902 0.2214 332.819 < 2e-16 ***
## as.factor(microsite)shrub.ambient 0.2471 0.4179 0.591 0.5544
## as.factor(microsite)shrub.soil 3.3369 0.4760 7.011 2.44e-12 ***
## as.factor(microsite)square -1.0049 0.4414 -2.277 0.0228 *
## as.factor(microsite)triangle -3.1949 0.5171 -6.179 6.56e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 483.3702)
##
## Null deviance: 10668889 on 21959 degrees of freedom
## Residual deviance: 10612393 on 21955 degrees of freedom
## AIC: 198057
##
## Number of Fisher Scoring iterations: 2
tab_model(lm.temp)
|
Â
|
temp
|
|
Predictors
|
Estimates
|
CI
|
p
|
|
(Intercept)
|
73.69
|
73.26 – 74.12
|
<0.001
|
|
microsite [shrub.ambient]
|
0.25
|
-0.57 – 1.07
|
0.554
|
|
microsite [shrub.soil]
|
3.34
|
2.40 – 4.27
|
<0.001
|
|
microsite [square]
|
-1.00
|
-1.87 – -0.14
|
0.023
|
|
microsite [triangle]
|
-3.19
|
-4.21 – -2.18
|
<0.001
|
|
Observations
|
21960
|
|
R2 Nagelkerke
|
0.924
|
library (emmeans)
## Welcome to emmeans.
## NOTE -- Important change from versions <= 1.41:
## Indicator predictors are now treated as 2-level factors by default.
## To revert to old behavior, use emm_options(cov.keep = character(0))
anova(lm.temp, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: temp
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 21959 10668889
## as.factor(microsite) 4 56496 21955 10612393 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans1 <- emmeans(lm.temp, pairwise~microsite)
lm.intensity <- glm(intensity~as.factor(microsite), data=shelter.shrub.open, family = "quasipoisson")
summary (lm.intensity)
##
## Call:
## glm(formula = intensity ~ as.factor(microsite), family = "quasipoisson",
## data = shelter.shrub.open)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -81.516 -56.469 -33.933 2.815 293.607
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.11118 0.01800 450.56 <2e-16 ***
## as.factor(microsite)shrub.ambient -0.71605 0.04520 -15.84 <2e-16 ***
## as.factor(microsite)shrub.soil -4.58886 0.41459 -11.07 <2e-16 ***
## as.factor(microsite)square -0.68732 0.04727 -14.54 <2e-16 ***
## as.factor(microsite)triangle -0.58210 0.05431 -10.72 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 6402.7)
##
## Null deviance: 62739687 on 12629 degrees of freedom
## Residual deviance: 54663384 on 12625 degrees of freedom
## (9330 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 6
tab_model(lm.intensity)
|
Â
|
intensity
|
|
Predictors
|
Incidence Rate Ratios
|
CI
|
p
|
|
(Intercept)
|
3331.50
|
3215.33 – 3450.44
|
<0.001
|
|
microsite [shrub.ambient]
|
0.49
|
0.45 – 0.53
|
<0.001
|
|
microsite [shrub.soil]
|
0.01
|
0.00 – 0.02
|
<0.001
|
|
microsite [square]
|
0.50
|
0.46 – 0.55
|
<0.001
|
|
microsite [triangle]
|
0.56
|
0.50 – 0.62
|
<0.001
|
|
Observations
|
12630
|
|
R2 Nagelkerke
|
1.000
|
anova(lm.intensity, test = "Chisq")
## Analysis of Deviance Table
##
## Model: quasipoisson, link: log
##
## Response: intensity
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 12629 62739687
## as.factor(microsite) 4 8076303 12625 54663384 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(lm.temp, pairwise~microsite)
## $emmeans
## microsite emmean SE df asymp.LCL asymp.UCL
## open 73.7 0.221 Inf 73.3 74.1
## shrub.ambient 73.9 0.354 Inf 73.2 74.6
## shrub.soil 77.0 0.421 Inf 76.2 77.9
## square 72.7 0.382 Inf 71.9 73.4
## triangle 70.5 0.467 Inf 69.6 71.4
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## open - shrub.ambient -0.247 0.418 Inf -0.591 0.9764
## open - shrub.soil -3.337 0.476 Inf -7.011 <.0001
## open - square 1.005 0.441 Inf 2.277 0.1525
## open - triangle 3.195 0.517 Inf 6.179 <.0001
## shrub.ambient - shrub.soil -3.090 0.551 Inf -5.612 <.0001
## shrub.ambient - square 1.252 0.521 Inf 2.403 0.1144
## shrub.ambient - triangle 3.442 0.586 Inf 5.869 <.0001
## shrub.soil - square 4.342 0.569 Inf 7.636 <.0001
## shrub.soil - triangle 6.532 0.629 Inf 10.382 <.0001
## square - triangle 2.190 0.603 Inf 3.629 0.0026
##
## P value adjustment: tukey method for comparing a family of 5 estimates
emmeans (lm.intensity, pairwise~microsite)
## $emmeans
## microsite emmean SE df asymp.LCL asymp.UCL
## open 8.111 0.01800 Inf 8.076 8.146
## shrub.ambient 7.395 0.04146 Inf 7.314 7.476
## shrub.soil 3.522 0.41420 Inf 2.711 4.334
## square 7.424 0.04371 Inf 7.338 7.510
## triangle 7.529 0.05124 Inf 7.429 7.630
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## open - shrub.ambient 0.7161 0.0452 Inf 15.842 <.0001
## open - shrub.soil 4.5889 0.4146 Inf 11.068 <.0001
## open - square 0.6873 0.0473 Inf 14.539 <.0001
## open - triangle 0.5821 0.0543 Inf 10.718 <.0001
## shrub.ambient - shrub.soil 3.8728 0.4163 Inf 9.304 <.0001
## shrub.ambient - square -0.0287 0.0602 Inf -0.477 0.9895
## shrub.ambient - triangle -0.1340 0.0659 Inf -2.032 0.2505
## shrub.soil - square -3.9015 0.4165 Inf -9.367 <.0001
## shrub.soil - triangle -4.0068 0.4174 Inf -9.600 <.0001
## square - triangle -0.1052 0.0674 Inf -1.562 0.5218
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates
#let's look at shape and blockage more closely
lm.temp.cover <- glm(temp~as.factor(microsite)* as.factor(cover.type), data = shelter.shrub.open, family="gaussian")
summary(lm.temp.cover)
##
## Call:
## glm(formula = temp ~ as.factor(microsite) * as.factor(cover.type),
## family = "gaussian", data = shelter.shrub.open)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -40.064 -16.940 -5.098 15.257 84.418
##
## Coefficients: (4 not defined because of singularities)
## Estimate Std. Error
## (Intercept) 73.6902 0.2168
## as.factor(microsite)square 30.1200 21.5514
## as.factor(microsite)triangle 26.5928 21.5320
## as.factor(cover.type)15 -26.4645 21.5505
## as.factor(cover.type)50 -30.1952 21.5505
## as.factor(cover.type)90 -31.2449 21.5405
## as.factor(microsite)square:as.factor(cover.type)15 -2.4960 1.4509
## as.factor(microsite)triangle:as.factor(cover.type)15 NA NA
## as.factor(microsite)square:as.factor(cover.type)50 -2.9435 1.4487
## as.factor(microsite)triangle:as.factor(cover.type)50 NA NA
## as.factor(microsite)square:as.factor(cover.type)90 NA NA
## as.factor(microsite)triangle:as.factor(cover.type)90 NA NA
## t value Pr(>|t|)
## (Intercept) 339.849 <2e-16 ***
## as.factor(microsite)square 1.398 0.1623
## as.factor(microsite)triangle 1.235 0.2168
## as.factor(cover.type)15 -1.228 0.2195
## as.factor(cover.type)50 -1.401 0.1612
## as.factor(cover.type)90 -1.451 0.1469
## as.factor(microsite)square:as.factor(cover.type)15 -1.720 0.0854 .
## as.factor(microsite)triangle:as.factor(cover.type)15 NA NA
## as.factor(microsite)square:as.factor(cover.type)50 -2.032 0.0422 *
## as.factor(microsite)triangle:as.factor(cover.type)50 NA NA
## as.factor(microsite)square:as.factor(cover.type)90 NA NA
## as.factor(microsite)triangle:as.factor(cover.type)90 NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 463.5796)
##
## Null deviance: 7168340 on 15388 degrees of freedom
## Residual deviance: 7130318 on 15381 degrees of freedom
## (6571 observations deleted due to missingness)
## AIC: 138155
##
## Number of Fisher Scoring iterations: 2
anova (lm.temp.cover, test="Chisq")#all variables are significant
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: temp
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df
## NULL 15388
## as.factor(microsite) 2 18912.1 15386
## as.factor(cover.type) 3 16739.0 15383
## as.factor(microsite):as.factor(cover.type) 2 2370.6 15381
## Resid. Dev Pr(>Chi)
## NULL 7168340
## as.factor(microsite) 7149427 1.384e-09 ***
## as.factor(cover.type) 7132688 7.104e-08 ***
## as.factor(microsite):as.factor(cover.type) 7130318 0.07755 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(lm.temp.cover, pairwise~microsite|cover.type)# square-triangle different at 15% and 90%
## $emmeans
## cover.type = 0:
## microsite emmean SE df asymp.LCL asymp.UCL
## open 73.7 0.217 Inf 73.3 74.1
## square nonEst NA NA NA NA
## triangle 100.3 21.531 Inf 58.1 142.5
##
## cover.type = 15:
## microsite emmean SE df asymp.LCL asymp.UCL
## open nonEst NA NA NA NA
## square 74.8 0.649 Inf 73.6 76.1
## triangle 73.8 0.920 Inf 72.0 75.6
##
## cover.type = 50:
## microsite emmean SE df asymp.LCL asymp.UCL
## open nonEst NA NA NA NA
## square 70.7 0.644 Inf 69.4 71.9
## triangle 70.1 0.920 Inf 68.3 71.9
##
## cover.type = 90:
## microsite emmean SE df asymp.LCL asymp.UCL
## open nonEst NA NA NA NA
## square 72.6 0.650 Inf 71.3 73.8
## triangle 69.0 0.644 Inf 67.8 70.3
##
## Confidence level used: 0.95
##
## $contrasts
## cover.type = 0:
## contrast estimate SE df z.ratio p.value
## open - square nonEst NA NA NA NA
## open - triangle -26.593 21.532 Inf -1.235 0.4324
## square - triangle nonEst NA NA NA NA
##
## cover.type = 15:
## contrast estimate SE df z.ratio p.value
## open - square nonEst NA NA NA NA
## open - triangle nonEst NA NA NA NA
## square - triangle 1.031 1.126 Inf 0.916 0.6301
##
## cover.type = 50:
## contrast estimate SE df z.ratio p.value
## open - square nonEst NA NA NA NA
## open - triangle nonEst NA NA NA NA
## square - triangle 0.584 1.123 Inf 0.520 0.8616
##
## cover.type = 90:
## contrast estimate SE df z.ratio p.value
## open - square nonEst NA NA NA NA
## open - triangle nonEst NA NA NA NA
## square - triangle 3.527 0.915 Inf 3.853 0.0003
##
## P value adjustment: tukey method for comparing a family of 3 estimates
#calculating Relative Interaction Index (RII)
#turn table into wide format
data.wide.temp <- reshape (shelter.shrub.open, timevar = "microsite", v.names = "temp", direction = "wide", idvar="date")
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar,
## varying = varying, : some constant variables (month,day,hour.
## 24,hour,lat,lon,rep,intensity,pendant.id,timeblock,soil.moisture,cover.type)
## are really varying
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=triangle: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=open: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=square: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.soil: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.ambient: first taken
data.wide.intensity<- reshape(shelter.shrub.open, timevar = "microsite", v.names = "intensity", direction= "wide", idvar = "date")
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar,
## varying = varying, : some constant variables (month,day,hour.
## 24,hour,lat,lon,rep,temp,pendant.id,timeblock,soil.moisture,cover.type) are
## really varying
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=triangle: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=open: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=square: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.soil: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.ambient: first taken
#rii for all temperature data
rii.temp.shrub<- data.wide.temp %>%
mutate(rii_calc_shrub = (temp.shrub.ambient-temp.open)/(temp.shrub.ambient + temp.open))
rii.temp.triangle<- data.wide.temp %>%
mutate(rii_calc_triangle = (temp.triangle-temp.open)/(temp.triangle + temp.open))
rii.temp.square<- data.wide.temp %>%
mutate(rii_calc_square = (temp.square-temp.open)/(temp.square+ temp.open))
#rii for all intensity data
rii.intensity.shrub<- data.wide.intensity %>%
mutate(rii_calc_shrub = (intensity.shrub.ambient-intensity.open)/(intensity.shrub.ambient+ intensity.open)) #it doesn't work as well as temperature
x <- select(rii.temp.shrub, rii_calc_shrub)
y <- select(rii.temp.triangle, rii_calc_triangle)
z <- select(rii.temp.square, rii_calc_square)
rii.final.temp<-cbind(x, y, z)
rii.csv<- read.csv("C:/Users/Nargol Ghazian/Desktop/Animal-Behaviour-and-Climate-project/CH3/RII.microsite.csv")
ggplot(rii.csv, aes(Microsite, rii)) +
geom_boxplot()+ ylab("Relative Interaction Index (RII)")+ geom_hline(yintercept=0, linetype="dashed", color = "red")+ theme_classic()

ggplot(rii.csv, aes(rii, fill = Microsite)) +
geom_histogram() +
geom_vline(xintercept = 0, col = 2, lty = 2) +
labs(x = "Relative Interaction Index (RII)", y = "Frequency", fill = "") +
scale_fill_brewer(palette = "Paired")+ theme_classic()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggline (rii.csv, x="Microsite", y= "rii", add = "mean_se", ylab = "Relative Interaction Index (RII)", xlab = "Microsite")

#GLM for RII
summary (rii.csv)
## Microsite rii
## shrub.ambient:14 Min. :-0.0409504
## square :24 1st Qu.:-0.0088791
## triangle :24 Median : 0.0006209
## Mean : 0.0025064
## 3rd Qu.: 0.0110682
## Max. : 0.0709935
lm.rii <- glm(rii~as.factor(Microsite), data = rii.csv, family="gaussian")
summary(lm.rii)
##
## Call:
## glm(formula = rii ~ as.factor(Microsite), family = "gaussian",
## data = rii.csv)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.044037 -0.009890 -0.002490 0.007168 0.067907
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.003086 0.005347 0.577 0.566
## as.factor(Microsite)square -0.006171 0.006728 -0.917 0.363
## as.factor(Microsite)triangle 0.004673 0.006728 0.695 0.490
##
## (Dispersion parameter for gaussian family taken to be 0.0004002093)
##
## Null deviance: 0.025029 on 61 degrees of freedom
## Residual deviance: 0.023612 on 59 degrees of freedom
## AIC: -304.19
##
## Number of Fisher Scoring iterations: 2
anova (lm.rii, test="Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: rii
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 61 0.025029
## as.factor(Microsite) 2 0.0014171 59 0.023612 0.1703
emmeans(lm.rii, pairwise~Microsite)#no significant difference
## $emmeans
## Microsite emmean SE df asymp.LCL asymp.UCL
## shrub.ambient 0.00309 0.00535 Inf -0.007393 0.01357
## square -0.00308 0.00408 Inf -0.011088 0.00492
## triangle 0.00776 0.00408 Inf -0.000245 0.01576
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## shrub.ambient - square 0.00617 0.00673 Inf 0.917 0.6294
## shrub.ambient - triangle -0.00467 0.00673 Inf -0.695 0.7667
## square - triangle -0.01084 0.00578 Inf -1.878 0.1451
##
## P value adjustment: tukey method for comparing a family of 3 estimates
library(sjPlot)
tab_model(lm.rii)
|
Â
|
rii
|
|
Predictors
|
Estimates
|
CI
|
p
|
|
(Intercept)
|
0.00
|
-0.01 – 0.01
|
0.566
|
|
Microsite [square]
|
-0.01
|
-0.02 – 0.01
|
0.363
|
|
Microsite [triangle]
|
0.00
|
-0.01 – 0.02
|
0.490
|
|
Observations
|
62
|
|
R2 Nagelkerke
|
0.057
|
ggplot (shelter.shrub.open, aes(as.factor(day), temp)) + geom_boxplot() + xlab("Day") + ylab ("Temperature (°F)")+ theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))+facet_grid("microsite")+ stat_summary(fun.y=mean, colour="red", geom="point", shape=18, size=1,show_guide = FALSE)
## Warning: `show_guide` has been deprecated. Please use `show.legend`
## instead.

ggplot (shelter.shrub.open, aes(as.factor(day), intensity)) + geom_boxplot() + xlab("Day") + ylab ("Solar Radiation (lum/ft²)")+ theme_classic()+ theme(axis.text=element_text(size=12))+ theme(axis.text=element_text(size=12))+facet_grid("microsite")+ stat_summary(fun.y=mean, colour="red", geom="point", shape=18, size=1,show_guide = FALSE)
## Warning: `show_guide` has been deprecated. Please use `show.legend`
## instead.
## Warning: Removed 9330 rows containing non-finite values (stat_boxplot).
## Warning: Removed 9330 rows containing non-finite values (stat_summary).

#let's create a heatmap
ggplot(shelter.shrub.open, aes(as.factor(day), as.factor(microsite), fill=temp))+ geom_tile()+ xlab("Day") + ylab ("Microsite")+theme(axis.text.x = element_text(angle = 90))+ scale_fill_distiller (palette='Blues')+ theme_classic()+ labs(fill = "Temperature (°F)")

#geom_jitter plot for max and min of intensity for each microsite
ggplot(shelter.shrub.open, aes(x = as.factor(microsite), y = intensity)) +
scale_y_log10() +
geom_jitter(position = position_jitter(width = 0.1, height = 0), alpha = 1/4) +
stat_summary(fun.y = min, colour = "blue", geom = "point", size = 5) +
stat_summary(fun.y = max, colour = "red", geom = "point", size = 5)+ xlab("")+ ylab("Solar Radiation (lum/ft²)")+theme_classic()
## Warning: Removed 9330 rows containing non-finite values (stat_summary).
## Warning: Removed 9330 rows containing non-finite values (stat_summary).
## Warning: Removed 9330 rows containing missing values (geom_point).

#faceted by microsite
ggplot(shelter.shrub.open, aes(x = as.factor(day), y = intensity)) +
scale_y_log10() +
geom_jitter(position = position_jitter(width = 0.1, height = 0), alpha = 1/4) +
stat_summary(fun.y = min, colour = "blue", geom = "point", size = 5) +
stat_summary(fun.y = max, colour = "red", geom = "point", size = 5)+ xlab("Day")+ ylab("Light Intensity (lum/ft²)")+theme_classic()+ facet_grid(~microsite)+ coord_flip()+ theme(axis.text.x = element_text(angle = 90))
## Warning: Removed 9330 rows containing non-finite values (stat_summary).
## Warning: Removed 9330 rows containing non-finite values (stat_summary).
## Warning: Removed 9330 rows containing missing values (geom_point).

#MACRO-MICRO CLIMATE CONTRAST
##test between macro-climate and micro-climate
library (ggplot2)
library(sjPlot)
ggplot(macro.micro.contrast, aes((microsite), temp, fill=microsite)) + geom_boxplot() + xlab("Site") + ylab ("Temperature (°F)")+ theme_classic()+ theme(axis.text=element_text(size=12))+theme(axis.text.x = element_text(angle = 90))+ theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())+ labs(fill = "Site")+ stat_summary(fun.y=mean, colour="black", geom="point", shape=18, size=3,show_guide = FALSE)
## Warning: `show_guide` has been deprecated. Please use `show.legend`
## instead.

#weather station data are significantly differnt from shrub, open, and square.
ggline(macro.micro.contrast, x = "day", y = "temp", color = "microsite",
add = "mean_ci", shape = "microsite", xlab = "Day", ylab = "Temperature (°F)", legend.title= "Site", legend="right")

ggline(macro.micro.contrast, x = "day", y = "intensity", color = "microsite",
add = "mean_ci", shape = "microsite", xlab = "Day", ylab = "Solar Rdation (lum/ft²)", legend.title= "Site", legend="right")
## Warning: Removed 9330 rows containing non-finite values (stat_summary).

lm.site <- glm(temp~as.factor(microsite), data = macro.micro.contrast, family="gaussian")
summary (lm.site)
##
## Call:
## glm(formula = temp ~ as.factor(microsite), family = "gaussian",
## data = macro.micro.contrast)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -39.543 -16.315 -5.344 13.475 110.218
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 73.6902 0.2193 335.970 < 2e-16
## as.factor(microsite)shrub.ambient 0.2471 0.4140 0.597 0.5506
## as.factor(microsite)shrub.soil 3.3369 0.4715 7.077 1.51e-12
## as.factor(microsite)square -1.0049 0.4373 -2.298 0.0216
## as.factor(microsite)triangle -3.1949 0.5122 -6.238 4.52e-10
## as.factor(microsite)weather.station -5.4202 0.8990 -6.029 1.68e-09
##
## (Intercept) ***
## as.factor(microsite)shrub.ambient
## as.factor(microsite)shrub.soil ***
## as.factor(microsite)square *
## as.factor(microsite)triangle ***
## as.factor(microsite)weather.station ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 474.3471)
##
## Null deviance: 10784021 on 22583 degrees of freedom
## Residual deviance: 10709809 on 22578 degrees of freedom
## AIC: 203260
##
## Number of Fisher Scoring iterations: 2
tab_model(lm.site)
|
Â
|
temp
|
|
Predictors
|
Estimates
|
CI
|
p
|
|
(Intercept)
|
73.69
|
73.26 – 74.12
|
<0.001
|
|
microsite [shrub.ambient]
|
0.25
|
-0.56 – 1.06
|
0.551
|
|
microsite [shrub.soil]
|
3.34
|
2.41 – 4.26
|
<0.001
|
|
microsite [square]
|
-1.00
|
-1.86 – -0.15
|
0.022
|
|
microsite [triangle]
|
-3.19
|
-4.20 – -2.19
|
<0.001
|
microsite [weather.station]
|
-5.42
|
-7.18 – -3.66
|
<0.001
|
|
Observations
|
22584
|
|
R2 Nagelkerke
|
0.963
|
library (emmeans)
anova(lm.site, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: temp
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 22583 10784021
## as.factor(microsite) 5 74212 22578 10709809 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(lm.site, pairwise~microsite)
## $emmeans
## microsite emmean SE df asymp.LCL asymp.UCL
## open 73.7 0.219 Inf 73.3 74.1
## shrub.ambient 73.9 0.351 Inf 73.2 74.6
## shrub.soil 77.0 0.417 Inf 76.2 77.8
## square 72.7 0.378 Inf 71.9 73.4
## triangle 70.5 0.463 Inf 69.6 71.4
## weather.station 68.3 0.872 Inf 66.6 70.0
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## open - shrub.ambient -0.247 0.414 Inf -0.597 0.9913
## open - shrub.soil -3.337 0.471 Inf -7.077 <.0001
## open - square 1.005 0.437 Inf 2.298 0.1947
## open - triangle 3.195 0.512 Inf 6.238 <.0001
## open - weather.station 5.420 0.899 Inf 6.029 <.0001
## shrub.ambient - shrub.soil -3.090 0.545 Inf -5.665 <.0001
## shrub.ambient - square 1.252 0.516 Inf 2.426 0.1473
## shrub.ambient - triangle 3.442 0.581 Inf 5.925 <.0001
## shrub.ambient - weather.station 5.667 0.940 Inf 6.030 <.0001
## shrub.soil - square 4.342 0.563 Inf 7.708 <.0001
## shrub.soil - triangle 6.532 0.623 Inf 10.480 <.0001
## shrub.soil - weather.station 8.757 0.967 Inf 9.059 <.0001
## square - triangle 2.190 0.598 Inf 3.664 0.0034
## square - weather.station 4.415 0.950 Inf 4.646 <.0001
## triangle - weather.station 2.225 0.987 Inf 2.254 0.2131
##
## P value adjustment: tukey method for comparing a family of 6 estimates
lm.site.sun <- glm(intensity~as.factor(microsite), data = macro.micro.contrast, family="gaussian")
summary (lm.site.sun)
##
## Call:
## glm(formula = intensity ~ as.factor(microsite), family = "gaussian",
## data = macro.micro.contrast)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -16478 -2436 -1303 124 47091
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3331.50 80.51 41.382 < 2e-16
## as.factor(microsite)shrub.ambient -1703.47 152.58 -11.165 < 2e-16
## as.factor(microsite)shrub.soil -3297.64 203.37 -16.215 < 2e-16
## as.factor(microsite)square -1656.02 160.31 -10.330 < 2e-16
## as.factor(microsite)triangle -1470.12 189.26 -7.768 8.58e-15
## as.factor(microsite)weather.station 13146.49 260.91 50.387 < 2e-16
##
## (Intercept) ***
## as.factor(microsite)shrub.ambient ***
## as.factor(microsite)shrub.soil ***
## as.factor(microsite)square ***
## as.factor(microsite)triangle ***
## as.factor(microsite)weather.station ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 38433720)
##
## Null deviance: 6.4238e+11 on 13253 degrees of freedom
## Residual deviance: 5.0917e+11 on 13248 degrees of freedom
## (9330 observations deleted due to missingness)
## AIC: 269095
##
## Number of Fisher Scoring iterations: 2
tab_model(lm.site.sun)
|
Â
|
intensity
|
|
Predictors
|
Estimates
|
CI
|
p
|
|
(Intercept)
|
3331.50
|
3173.71 – 3489.29
|
<0.001
|
|
microsite [shrub.ambient]
|
-1703.47
|
-2002.52 – -1404.43
|
<0.001
|
|
microsite [shrub.soil]
|
-3297.64
|
-3696.23 – -2899.05
|
<0.001
|
|
microsite [square]
|
-1656.02
|
-1970.22 – -1341.83
|
<0.001
|
|
microsite [triangle]
|
-1470.12
|
-1841.06 – -1099.17
|
<0.001
|
microsite [weather.station]
|
13146.49
|
12635.12 – 13657.86
|
<0.001
|
|
Observations
|
13254
|
|
R2 Nagelkerke
|
1.000
|
library (emmeans)
anova(lm.site.sun, test = "Chisq")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: intensity
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 13253 6.4238e+11
## as.factor(microsite) 5 1.3321e+11 13248 5.0917e+11 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(lm.site.sun, pairwise~microsite)
## $emmeans
## microsite emmean SE df asymp.LCL asymp.UCL
## open 3331.5 80.5 Inf 3174 3489
## shrub.ambient 1628.0 129.6 Inf 1374 1882
## shrub.soil 33.9 186.8 Inf -332 400
## square 1675.5 138.6 Inf 1404 1947
## triangle 1861.4 171.3 Inf 1526 2197
## weather.station 16478.0 248.2 Inf 15992 16964
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## open - shrub.ambient 1703.5 153 Inf 11.165 <.0001
## open - shrub.soil 3297.6 203 Inf 16.215 <.0001
## open - square 1656.0 160 Inf 10.330 <.0001
## open - triangle 1470.1 189 Inf 7.768 <.0001
## open - weather.station -13146.5 261 Inf -50.387 <.0001
## shrub.ambient - shrub.soil 1594.2 227 Inf 7.013 <.0001
## shrub.ambient - square -47.5 190 Inf -0.250 0.9999
## shrub.ambient - triangle -233.4 215 Inf -1.086 0.8870
## shrub.ambient - weather.station -14850.0 280 Inf -53.039 <.0001
## shrub.soil - square -1641.6 233 Inf -7.058 <.0001
## shrub.soil - triangle -1827.5 253 Inf -7.212 <.0001
## shrub.soil - weather.station -16444.1 311 Inf -52.944 <.0001
## square - triangle -185.9 220 Inf -0.844 0.9593
## square - weather.station -14802.5 284 Inf -52.072 <.0001
## triangle - weather.station -14616.6 302 Inf -48.472 <.0001
##
## P value adjustment: tukey method for comparing a family of 6 estimates
ggplot(macro.micro.contrast, aes(temp, fill = microsite)) +
geom_histogram(binwidth = 5) +
scale_fill_brewer(palette = "Set1")+ labs(fill = "", x = "Temperature (°F)", y = "Frequency")+theme_classic()+ theme(axis.text=element_text(size=12))+ stat_central_tendency(aes(color = microsite), type = "mean", color="green", linetype = 2)+ stat_central_tendency(type = "median", color = "blue", linetype = 2)+ facet_grid(~microsite)+theme(legend.position = "none")

ggplot(macro.micro.contrast, aes((day), temp, color=microsite)) + geom_smooth()+ xlab("Day") + ylab ("Temperature (°F)")+ theme_classic()+ theme(axis.text=element_text(size=12))+stat_summary(fun.y=max, geom="point", size=2, aes(shape = microsite))+ facet_grid("microsite")+ theme(legend.position = "none")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

ggplot(macro.micro.contrast, aes((day), intensity, color=microsite)) + geom_smooth()+ xlab("Day") + ylab ("Solar Rdation (lum/ft²)")+ theme_classic()+ theme(axis.text=element_text(size=12))+stat_summary(fun.y=max, geom="point", size=2, aes(shape = microsite))+ facet_grid("microsite")+ theme(legend.position = "none")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 9330 rows containing non-finite values (stat_smooth).
## Warning: Removed 9330 rows containing non-finite values (stat_summary).

library(car)
## Warning: package 'car' was built under R version 3.6.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 3.6.3
## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
leveneTest(temp ~ microsite, macro.micro.contrast)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 60.096 < 2.2e-16 ***
## 22578
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
leveneTest(intensity ~microsite, macro.micro.contrast)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 815.31 < 2.2e-16 ***
## 13248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#variation in sunlight and intensity are significant
wide.clim <- reshape(macro.micro.contrast, v.names="temp", timevar="microsite", idvar=c("temp", "microsite"),
direction="wide")
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar,
## varying = varying, : some constant variables (date,month,day,hour.
## 24,hour,lat,lon,rep,intensity,pendant.id,timeblock,soil.moisture,cover.type)
## are really varying
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=triangle: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=open: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=square: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.soil: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.ambient: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=weather.station: first taken
summary(wide.clim)
## site sensor date year
## panoche:2671 hobo.pendant:2429 05/27/19: 168 Min. :2019
## satellite : 242 05/24/19: 133 1st Qu.:2019
## 05/21/19: 130 Median :2019
## 05/20/19: 126 Mean :2019
## 05/25/19: 126 3rd Qu.:2019
## 05/29/19: 123 Max. :2019
## (Other) :1865
## month day hour.24 hour
## Min. :5.000 Min. : 1.0 10:46:02: 22 : 242
## 1st Qu.:5.000 1st Qu.: 7.0 8:22:41 : 22 10:46:02 AM: 22
## Median :5.000 Median :21.0 1:38:16 : 21 8:22:41 AM : 22
## Mean :5.364 Mean :17.6 10:22:04: 21 1:38:16 AM : 21
## 3rd Qu.:6.000 3rd Qu.:26.0 10:38:01: 20 10:22:04 AM: 21
## Max. :6.000 Max. :31.0 12:38:43: 20 1:24:11 PM : 20
## (Other) :2545 (Other) :2323
## lat lon rep intensity
## Min. :36.39 Min. :-120.9 Min. : 1.000 Min. : 0
## 1st Qu.:36.69 1st Qu.:-120.8 1st Qu.: 1.000 1st Qu.: 96
## Median :36.69 Median :-120.8 Median : 2.000 Median : 768
## Mean :36.69 Mean :-120.8 Mean : 2.404 Mean : 5697
## 3rd Qu.:36.70 3rd Qu.:-120.8 3rd Qu.: 3.000 3rd Qu.: 4864
## Max. :36.70 Max. :-120.8 Max. :12.000 Max. :63505
## NA's :242 NA's :242 NA's :242 NA's :727
## pendant.id timeblock soil.moisture cover.type
## Min. : 1156629 :1024 Min. : 1.00 Min. : 0.00
## 1st Qu.: 1174238 morning : 714 1st Qu.: 23.00 1st Qu.: 0.00
## Median :20567072 evening : 563 Median : 29.40 Median :15.00
## Mean :13343324 afternoon: 128 Mean : 82.71 Mean :21.51
## 3rd Qu.:20567073 0 : 50 3rd Qu.: 31.80 3rd Qu.:15.00
## Max. :20619223 0.0001 : 3 Max. :998.00 Max. :90.00
## NA's :242 (Other) : 189 NA's :1074 NA's :1266
## temp.triangle temp.open temp.square temp.shrub.soil
## Min. : 36.53 Min. : 34.20 Min. : 34.79 Min. : 45.10
## 1st Qu.: 56.84 1st Qu.: 58.22 1st Qu.: 58.56 1st Qu.: 66.62
## Median : 74.62 Median : 80.21 Median : 79.15 Median : 87.55
## Mean : 75.37 Mean : 81.84 Mean : 80.63 Mean : 91.82
## 3rd Qu.: 93.21 3rd Qu.:104.23 3rd Qu.:101.38 3rd Qu.:113.47
## Max. :124.08 Max. :156.98 Max. :156.98 Max. :187.25
## NA's :2256 NA's :2160 NA's :2192 NA's :2192
## temp.shrub.ambient temp.weather.station
## Min. : 34.39 Min. : 43.88
## 1st Qu.: 60.54 1st Qu.: 58.87
## Median : 84.11 Median : 69.89
## Mean : 87.14 Mean : 70.90
## 3rd Qu.:110.92 3rd Qu.: 82.72
## Max. :162.72 Max. :100.58
## NA's :2126 NA's :2429
wide.clim.sun <- reshape(macro.micro.contrast, v.names="intensity", timevar="microsite", idvar=c("temp", "microsite"), direction="wide")
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar,
## varying = varying, : some constant variables (date,month,day,hour.
## 24,hour,lat,lon,rep,pendant.id,timeblock,soil.moisture,cover.type) are
## really varying
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=triangle: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=open: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=square: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.soil: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=shrub.ambient: first taken
## Warning in reshapeWide(data, idvar = idvar, timevar = timevar, varying =
## varying, : multiple rows match for microsite=weather.station: first taken
summary (wide.clim.sun)
## site sensor date year
## panoche:2671 hobo.pendant:2429 05/27/19: 168 Min. :2019
## satellite : 242 05/24/19: 133 1st Qu.:2019
## 05/21/19: 130 Median :2019
## 05/20/19: 126 Mean :2019
## 05/25/19: 126 3rd Qu.:2019
## 05/29/19: 123 Max. :2019
## (Other) :1865
## month day hour.24 hour
## Min. :5.000 Min. : 1.0 10:46:02: 22 : 242
## 1st Qu.:5.000 1st Qu.: 7.0 8:22:41 : 22 10:46:02 AM: 22
## Median :5.000 Median :21.0 1:38:16 : 21 8:22:41 AM : 22
## Mean :5.364 Mean :17.6 10:22:04: 21 1:38:16 AM : 21
## 3rd Qu.:6.000 3rd Qu.:26.0 10:38:01: 20 10:22:04 AM: 21
## Max. :6.000 Max. :31.0 12:38:43: 20 1:24:11 PM : 20
## (Other) :2545 (Other) :2323
## lat lon rep temp
## Min. :36.39 Min. :-120.9 Min. : 1.000 Min. : 34.20
## 1st Qu.:36.69 1st Qu.:-120.8 1st Qu.: 1.000 1st Qu.: 59.88
## Median :36.69 Median :-120.8 Median : 2.000 Median : 79.50
## Mean :36.69 Mean :-120.8 Mean : 2.404 Mean : 82.50
## 3rd Qu.:36.70 3rd Qu.:-120.8 3rd Qu.: 3.000 3rd Qu.:101.69
## Max. :36.70 Max. :-120.8 Max. :12.000 Max. :187.25
## NA's :242 NA's :242 NA's :242
## pendant.id timeblock soil.moisture cover.type
## Min. : 1156629 :1024 Min. : 1.00 Min. : 0.00
## 1st Qu.: 1174238 morning : 714 1st Qu.: 23.00 1st Qu.: 0.00
## Median :20567072 evening : 563 Median : 29.40 Median :15.00
## Mean :13343324 afternoon: 128 Mean : 82.71 Mean :21.51
## 3rd Qu.:20567073 0 : 50 3rd Qu.: 31.80 3rd Qu.:15.00
## Max. :20619223 0.0001 : 3 Max. :998.00 Max. :90.00
## NA's :242 (Other) : 189 NA's :1074 NA's :1266
## intensity.triangle intensity.open intensity.square intensity.shrub.soil
## Min. : 2 Min. : 1 Min. : 1 Min. : 1.00
## 1st Qu.: 384 1st Qu.: 468 1st Qu.: 375 1st Qu.: 4.50
## Median : 1184 Median : 1752 Median : 1184 Median : 16.00
## Mean : 3938 Mean : 5631 Mean : 2376 Mean : 53.68
## 3rd Qu.: 5568 3rd Qu.: 9216 3rd Qu.: 3712 3rd Qu.: 29.00
## Max. :19456 Max. :27648 Max. :24576 Max. :1152.00
## NA's :2381 NA's :2299 NA's :2329 NA's :2380
## intensity.shrub.ambient intensity.weather.station
## Min. : 1 Min. : 0.0
## 1st Qu.: 248 1st Qu.: 334.1
## Median : 864 Median :17467.0
## Mean : 2476 Mean :24799.5
## 3rd Qu.: 3968 3rd Qu.:48678.7
## Max. :19456 Max. :63504.9
## NA's :2264 NA's :2429
#Create microsite map for appendix
library(ggmap)
## Warning: package 'ggmap' was built under R version 3.6.3
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
##
## Attaching package: 'ggmap'
## The following object is masked from 'package:magrittr':
##
## inset
register_google(key="AIzaSyBpfKtYrkYVS3LEJSjV1cIHeYrxJPsPX4U")
panoche1 <- get_map(location = c(lon = -120.7932, lat = 36.69363), zoom = 6, maptype = "satellite")
## Source : https://maps.googleapis.com/maps/api/staticmap?center=36.69363,-120.7932&zoom=6&size=640x640&scale=2&maptype=satellite&language=en-EN&key=xxx
panoche2 <- get_map(location = c(lon = -120.7932, lat = 36.69363), zoom = 12, maptype = "terrain")
## Source : https://maps.googleapis.com/maps/api/staticmap?center=36.69363,-120.7932&zoom=12&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx
lat_long <- read.csv("C:/Users/Nargol Ghazian/Desktop/Animal-Behaviour-and-Climate-project/CH3/Shrub_Shelter_Open lat long.csv")
na.omit(lat_long)
## Shelter.ID Lat Long Microsite
## 1 1 36.69363 -120.7932 Triangle
## 2 2 36.69364 -120.7933 Square
## 3 3 36.69355 -120.7931 Square
## 4 4 36.69349 -120.7932 Triangle
## 5 5 36.69349 -120.7931 Triangle
## 6 6 36.39342 -120.7931 Square
## 7 7 36.69394 -120.7930 Square
## 8 8 36.69397 -120.7929 Triangle
## 9 9 36.69401 -120.7928 Square
## 10 10 36.69400 -120.7930 Triangle
## 11 11 36.69405 -120.7930 Square
## 12 12 36.69408 -120.7930 Triangle
## 13 1 36.69534 -120.7970 Shrub
## 14 2 36.69532 -120.7972 Shrub
## 15 3 36.69534 -120.7972 Shrub
## 16 4 36.69606 -120.7967 Shrub
## 17 5 36.69593 -120.7967 Shrub
## 18 6 36.69596 -120.7968 Shrub
## 19 7 36.69589 -120.7980 Shrub
## 20 1 36.69529 -120.7970 Open
## 21 2 36.69533 -120.7971 Open
## 22 3 36.69529 -120.7923 Open
## 23 4 36.69596 -120.7966 Open
## 24 5 36.69587 -120.7967 Open
## 25 6 36.69601 -120.7967 Open
## 26 7 36.69590 -120.7981 Open
ggmap(panoche2) +
geom_point(data= lat_long, aes(x=Long, y=Lat, color = Microsite), alpha = 6/10, size =2, show.legend = TRUE) +
labs(title = "Map of Shelters and Shrubs", x = "Longitude", y = "Latitude")
## Warning: Removed 1 rows containing missing values (geom_point).

ggmap(panoche1) +
geom_point(data= lat_long, aes(x=Long, y=Lat, color = Microsite), alpha = 6/10, size =2, show.legend = TRUE) +
labs(title = "Map of Califronia", x = "Longitude", y = "Latitude")

#Very general arial view of microsites
#CONCLUSIONS
##1.Temperature and sunlight intensity are positively related.
##2.Shelters function similar to shrubs, but square emulates shrub slightly better.
##3.Triangle at 90% blockage is best at lowering temperature.
##4.The open experienced the most variation.
##5. Temperature from Weather station data were significantly lower than square, shrub, and the open- weather-station data is nor ecologically-relevant.